 
        %==========================================================
        % Simulate the economy
        % =============== ==========================================
        % This lists all possible combinations of individual income
        % values (IDIO_DIST), and the aggregate states (pair of individual income and ROV-income) associated
        % with every individual income value
        
        stable=(find(STABLE==1));
        
        if coaldev==1                       % for coalitional deviations, all the stable ones
            vss_set_sim=stable;
        else                               
            vss_set_sim=vss_set;
            
        end
        
            
            state=STATE(vss);
            if rov_average==1
                ScaleRov=vss;
            else
                ScaleRov=1;
            end
            
            S=S_KEEP(1:state,:,vss);
            S_rov=S_rov_KEEP(1:state/length(incomevals),:,vss);
            idio_dist=IDIO_DIST_ROV_KEEP(1:state/length(incomevals),1:length(incomevals),vss);
            phisol=PHI_KEEP(:,1:state,:,vss);
            
            if length(incomevals)==3
            IDIO_DIST=[ones(state/3,1)*[1 0 0] + idio_dist; ones(state/3,1)*[0 1 0] + idio_dist; ones(state/3,1)*[0 0 1]+idio_dist];
            elseif length(incomevals)==2
            IDIO_DIST=[ones(state/2,1)*[1 0] + idio_dist; ones(state/2,1)*[0 1] + idio_dist];
            end
            distribution=size(IDIO_DIST,1);
            
            
            var_dyagglog=zeros(Nvill*(vss+1),vss_high);var_dyvillagglog=zeros(Nvill*(vss+1),vss_high);vardaggylog=zeros(Nvill*(vss+1),vss_high);betayc=zeros(2,Nvill*(vss+1));vardcshare=zeros(Nvill*(vss+1),vss_high);vardcshare_dylow=zeros(Nvill*(vss+1),vss_high);vardcshare_dyhigh=zeros(Nvill*(vss+1),vss_high);varclog=zeros(Nvill*(vss+1),vss_high);vardclog=zeros(Nvill*(vss+1),vss_high);varylog=zeros(Nvill*(vss+1),vss_high);vardylog=zeros(Nvill*(vss+1),vss_high);varclog_cond=zeros(Nvill*(vss+1),vss_high);varylog_cond=zeros(Nvill*(vss+1),vss_high);
            varc_yhigh_cond=zeros(Nvill*(vss+1),vss_high);varc_ylow_cond=zeros(Nvill*(vss+1),vss_high);vardclog_cond=zeros(Nvill*(vss+1),vss_high);vardylog_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high);
            vardc_dypos_cond=zeros(Nvill*(vss+1),vss_high);
            vardc_dyneg_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high); varc_ylow=zeros(Nvill*(vss+1),vss_high);vardc_dypos=zeros(Nvill*(vss+1),vss_high); vardc_dyneg=zeros(Nvill*(vss+1),vss_high);
            
            betadcdy_agg=zeros(3,Nvill*(vss+1));   betadcdy=zeros(2,vss_high);  betadcdy_high=zeros(1,Nvill*(vss+1));betadcdy_low=zeros(1,Nvill*(vss+1));relcvar=zeros(Nvill*(vss+1),vss_high); reldcvar=zeros(Nvill*(vss+1),vss_high); reldc0var=zeros(Nvill*(vss+1),vss_high); relcvar_cond=zeros(Nvill*(vss+1),vss_high); reldcvar_cond=zeros(Nvill*(vss+1),vss_high);
            beta_cond=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_cond=zeros(4,Nvill*(vss+1),vss_high); beta_yhigh_cond=zeros(4,Nvill*(vss+1),vss_high) ;
            beta_agg=zeros(2,Nvill*(vss+1),vss_high); beta_yhigh_agg=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_agg=zeros(2,Nvill*(vss+1),vss_high);
            
            
            %=============================================
            % Start simulation
            % ============================================
            % number of insurance groups equals village size
            % vss_high+1 divided by size of ins group vss+1
            
            %stream = RandStream.getGlobalStream;
            rng(3)
            Npanel=Nvill*(vss+1);
            csim0=nan(Nvill,vss+1,T);
            income0=nan(Nvill,vss+1,T);
            income0ind=nan(Nvill,vss+1,T);
            phisim0=nan(Nvill,vss+1,T);
            Yvill=nan(1,T);
            Yagg=nan(Nvill,T);
            Ssim=nan(Nvill,vss+1,T);
            gammarfun=nan(Nvill,vss+1,T);
            phinew=nan(Nvill,vss+1,T);
            indSjj=nan(T,vss+1);
            
            disp('simulating')
            
            disp('Reminder: now simulating for coaldev, exact')
            disp([coaldev exact])
            disp('now simulating for parameters')
            disp([Qbeta,Qrho])
            disp('vss')
            disp(vss)
            
            
            tic
            
            %initial state
            t=1;
            % initial values of income and consumption (set =
            % income)
            if length(incomevals)==3
            rng(vss,'twister');    
            shockinit=randi(3,Nvill,vss+1);
            elseif length(incomevals)==2
            rng(vss,'twister');    
            shockinit=randi(2,Nvill,vss+1);
            end
            income0ind(:,:,1)=shockinit;
            income0(:,:,1)=incomevals(income0ind(:,:,1));
            csim0(:,:,1)=income0(:,:,1);
            
            
            Yagg(:,t)=sum(income0(:,:,t),2);
            
            
            IDIO_DIST_sim=[];
            
            while t<T
                t=t+1;
                %%%%simulate next period income taking into account
                %%%% persistence across time
                rng(t,'twister');
                shock=rand(Nvill,vss+1);
                for k=1:vss+1
                    if length(incomevals)==3
                    rr=(find(shock(:,k)<=cumQ(income0ind(:,k,t-1),1)));
                    income0ind(rr,k,t)=1;
                    rr=(find(cumQ(income0ind(:,k,t-1),1)<shock(:,k) & cumQ(income0ind(:,k,t-1),2)>=shock(:,k)));
                    income0ind(rr,k,t)=2;
                    rr=(find(cumQ(income0ind(:,k,t-1),2)<shock(:,k) & cumQ(income0ind(:,k,t-1),3)>=shock(:,k)));
                    income0ind(rr,k,t)=3;
                    income0(:,k,t)=incomevals(income0ind(:,k,t));
                    elseif length(incomevals)==2
                    rr=(find(shock(:,k)<=cumQ(income0ind(:,k,t-1),1)));
                    income0ind(rr,k,t)=1;
                    rr=(find(cumQ(income0ind(:,k,t-1),1)<shock(:,k) & cumQ(income0ind(:,k,t-1),2)>=shock(:,k)));
                    income0ind(rr,k,t)=2;
                    income0(:,k,t)=incomevals(income0ind(:,k,t));
                    end
                end
               
                
                Yagg(:,t)=sum(income0(:,:,t),2);
                gammarfun(:,:,t)=(Yagg(:,t-1)*ones(1,vss+1)-csim0(:,:,t-1))./csim0(:,:,t-1)/ScaleRov;
                
                
                for ivill=1:Nvill
                    if length(incomevals)==3
                    IDIO_DIST_sim=find(((IDIO_DIST(:,1,1)==(sum(income0ind(ivill,:,t)==1))).*(IDIO_DIST(:,2,1)==(sum(income0ind(ivill,:,t)==2))).*(IDIO_DIST(:,3,1)==(sum(income0ind(ivill,:,t)==3))))==1)';
                    elseif length(incomevals)==2
                    IDIO_DIST_sim=find(((IDIO_DIST(:,1,1)==(sum(income0ind(ivill,:,t)==1))).*(IDIO_DIST(:,2,1)==(sum(income0ind(ivill,:,t)==2))))==1)';
                    end
                    for k=1:vss+1
                        % choose for individual k the aggregate state
                        % that applies
                        indSjj(t,k)=IDIO_DIST_sim(S(IDIO_DIST_sim,1)==income0(ivill,k,t));
                        % check that aggregate and individual income0s are
                        % consistent with the states caculated above
                        if abs(Yagg(ivill,t)-(S(indSjj(t,k),1)+ScaleRov*S(indSjj(t,k),2)))>0.000001
                            stop
                        end
                      
                       
                    end
                    Ssim(ivill,:,t)=indSjj(t,:);
                end
                
                
                for k=1:vss+1
                    BLA=interp1(gammagrid(:,4),phisol(:,Ssim(:,k,t),1),gammarfun(:,k,t));
                    phinew(:,k,t)=diag(BLA);
                end
                
               
                gammaraux=((1+phinew(:,:,t))./(1+phinew(:,1,t)*ones(1,vss+1))).^(1/rho).*csim0(:,:,t-1)./(csim0(:,1,t-1)*ones(1,vss+1));
                
                
                csim0(:,:,t)=(gammaraux./(sum(gammaraux,2)*ones(1,vss+1))).*(Yagg(:,t)*ones(1,vss+1));
                %end
            end
            disp('simulating takes')
            tsim(Qbeta,Qrho)=toc
            
            %%%%save the income simulation 
           
            save('INCOMESIM.mat','income0','income0ind');
             
            tic
            %reshape
            csim0keep=csim0;
            income0keep=income0;
            AA=[];
            BB=[];
            CC=[];
            for t=1:T
                AA=[AA reshape(csim0(:,:,t)',Nvill*(vss+1),1)];
                BB=[BB reshape(income0(:,:,t)',Nvill*(vss+1),1)];
                CC=[CC reshape(phinew(:,:,t)',Nvill*(vss+1),1)];
            end
            csim0=AA;
            income0=BB;
            phinew0=CC;
            disp('reshaping takes')
            toc

            
           